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Abstract —  In  imaging  applications  based  on  RF  Tomography, 
coupling  between  transmitters  and  receivers  is  a  principal 
technical  challenge  to  be  overcome.  To  mitigate  the  coupling,  we 
propose  to  activate  different  transmitters  simultaneously,  and 
determine  their  respective  current  source  distribution  in  order 
to  create  electric  field  nulls  at  desired  points  (e.g.,  the  receivers). 
The  current  design  problem  must  take  into  consideration 
physical  and  system  constraints  of  the  radiators.  As  such,  the 
current  design  must  be  formulated  as  a  nonlinear  programming 
problem.  These  constraints  are  discussed  in  details,  and  our 
method  is  validated  using  FDTD  simulations. 

I.  Introduction 

Radio  frequency  (RF)  tomography  is  an  imaging  technique 
in  which  distributed  transmitters  emit  narrowband  low- 
frequency  signals  into  the  area  of  interest,  and  distributed 
receivers  collect  the  scattered  field  from  targets  [1],  [7]. 
These  spatially  diverse  measurements  are  then  used  to  form 
an  image  of  the  belowground  scene.  At  low  frequencies,  both 
the  transmit  (Tx)  and  receive  (Rx)  radiation  patterns  are  well 
represented  as  electrically  small  dipoles  /  loops,  which  afford 
a  very  accurate  linear  approximation  to  the  forward  scattering 
model  [1].  As  such,  a  plethora  of  inversion  procedures  can  be 
used  to  produce  images,  even  with  subwavelength  and  range- 
independent  resolution.  This  is  true  even  for  very  large 
areas/volumes  of  interest,  both  in  the  near  and  the  far  field 
regions.  Hitherto,  RF  Tomography  has  been  considered 
mainly  for  belowground  imaging  of  tunnels  and  underground 
facilities  (UGF)  [1]  (see  Fig.  1);  however,  this  technique 
might  also  be  applied  to  close-in  sensing  of  urban 
environments,  multistatic  radar  imaging,  foliage  penetration, 
building  penetration,  etc.  [11]. 

Unlike  wideband  imaging  techniques  (e.g.,  time  reversal, 
travel-time  tomography),  which  are  based  on  transmission  of 
short  pulses,  RF  Tomography  cannot  use  time  gating  to 
isolate  the  direct  path  between  Tx  and  Rx  antennas  [1]. 
Furthermore,  returns  from  scattering  bodies  that  are  not  of 
interest  (i.e.,  clutter)  can  mask  the  return  from  weaker 
scatterers  that  are  of  interest.  In  principle,  if  the  background 
field  could  be  measured  prior  to  the  insertion  of  the  objects  of 
interest,  then  the  background  clutter  could  be  subtracted  [7], 


[8].  However,  this  approach  is  not  practical  for  many 
applications,  particularly  for  UGF  detection. 

The  contribution  of  this  paper  is  a  technique  for  mitigating 
the  direct  path  and  the  strong  return  from  known  scatterers  in 
RF  tomography  applications.  This  technique  is  based  on 
multiple  transmitters  simultaneously  emitting  signals  with 
source  currents  designed  to  create  electric  field  nulls  at 
receiver  and  clutter  locations. 

The  paper  is  structured  as  follows.  In  Section  II  we 
introduce  a  mathematical  description  of  the  geometry  used  in 
RF  Tomography.  In  Section  III  we  derive  the  3D  forward 
scattering  model  for  the  case  of  simultaneously  emitting 
transmitters.  In  Section  IV,  we  discuss  the  source  current 
design  problem,  which  is  formulated  as  a  nonlinear 
optimization  program.  In  Section  V,  we  validate  our  approach 
using  an  FDTD  simulation.  Concluding  remarks  are  given  in 
Section  VI. 


Figure  1:  pictorial  representation  of  RF  Tomography  for 
close-in  sensing  of  underground  facilities. 

II.  Overview  of  RF  Tomography 

Consider  the  notional  3D  transmitter-receiver  geometry 
depicted  in  Fig.  2.  The  /2-th  transmitter  is  located  at  position 

T and  the  /22-th  receiver  is  located  at  position  Xhm  .  In  order  to 
simplify  the  analysis,  we  assume  that  each  transmitter  emits  a 
monochromatic  signal  with  frequency  f .  The  multi- 
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frequency  case  can  be  derived  by  invoking  the  principle  of 
superposition.  We  further  simplify  the  analysis  by  assuming 
the  air-earth  interface  to  be  flat:  in  principle,  our  approach 
could  be  extended  to  irregular  surfaces  as  well. 

The  air  semispace  is  modeled  as  free-space  medium,  while 
the  ground  semispace  is  modeled  as  an  homogeneous 

medium  with  background  relative  dielectric  permittivity  £D  , 
background  conductivity  Od  ,  and  magnetic  permeability 

//0 .  The  dielectric/conducting  targets  (e.g.,  tunnels)  are 

assumed  to  reside  in  the  investigation  domain  D.  Transmitters 
and  receivers  must  reside  outside  the  investigation  domain 
D .  A  scattering  body  located  at  position  r f  inside  D  is 
described  by  a  deviation  in  the  relative  dielectric  permittivity 

£r(r?)  and  the  conductivity  Unlike  in  other 

geophysical  applications,  the  detection  of  underground 
facilities  does  not  require  the  discrimination  between 
dielectric  and  conducting  bodies.  Hence,  we  can  for 
convenience  define  a  complex- valued  contrast  function : 

V{r')  =  £r(r')-£D+j{(j(Y')-(jD)l2nf£0,  (1) 

where  £{)  is  the  dielectric  permittivity  of  free  space,  and 
j2  =  —  1 .  A  scattering  body  at  position  r ?  can  be  detected 
by  estimating  V (rf)  from  measurements  of  the  scattered 
field,  and  comparing  its  magnitude  to  a  threshold.  In  the 
following  section,  we  relate  to  the  measurements  of 

the  scattering  field. 


Figure  2:  three  dimensional  geometry  for  the  model. 


III.  Forward  Model 

To  properly  address  the  coupling  problem  we  must  account 
for  the  3D  vector  nature  of  the  electric  field  in  the 
development  of  the  forward  scattering  model.  We  assume  that 
the  ft-th  transmitter  is  constituted  of  (up  to)  three  co-located 
but  orthogonal  dipole  antennas.  Oftentimes,  the  orientation  of 
the  transmitters  and  receivers  will  not  be  congruent  with  the 
intertial  coordinate  system.  To  account  for  the  actual 


orientation,  we  assume  the  first,  second  and  third  dipole  of 
the  72-th  transmitter  to  be  oriented  along  the  unit  vector 

^nXn->  I*  >  respectively,  where  !~n  £n  and|w  correspond  to 

the  inertial  coordinate  vectors  X ,  y  and  Z  rotated  by  OCn , 

[5n  and  yn ,  respectively.  This  operation  can  be  represented 

by  the  rotation  matrix  ,  which  is  given  by: 

cos  y  cos  B  -sin  v  cos  B  sin  B 

"  "  "  "  (2) 

Rn  =  sin  yn  cos  OCn  +cos  Jn  sin/?  sin  OCn  cos  yn  cos  OCn  -  sin  Jn  sin  /?  sin  C(n  -cos  A  sin  dr 

sin  y  sin  a  -  cos  y  sin  f3  cos  CCn  cos  y  sin  a  +  sin  y  sin  cos  a  cos  /?  cos  0£n 

Currents  flowing  in  each  orthogonal  radiating  dipole  of  the  n- 
th  transmitter  can  be  expressed  in  phasor  form  by  the  vector: 

d«  =(dl  +Jdi  )?»  +{dl„  +Jdi  )L  +{dl  +Jdl  )l»  (3) 

Clearly,  from  the  phasor  dn  the  current  to  be  fed  into  each 
orthogonal  radiating  dipole  can  be  computed  as: 

d  (0  =  Re(d  exp(-y^)) 

=  Re(dw)cos(^)  +  Im(dw)sin(^)  (4) 

=  dgn  (t)in+dCn  (t)L+d^n  (0  In 

dg  (t)  =  ^i^n  Sin(^  +  tan  1  (dg  !  dg  ))  (5) 

( t )  =  yj(d^  )  +(df  )  sin^^  +  tan-1  (d^  / dJc  (6) 

dg  (0  =  ^(dg  )  +(d^  )  sin {^cat  +  tan-1  / dJ^  ))  (7) 

where  CO  =  'll!  f  .  Clearly,  any  desired  phasor  form  dn  can 
be  synthesized  by  driving  the  dipoles  with  the  appropriate 
current  amplitude  and  phase.  The  phasor  dn,  which  is 

measured  with  respect  to  a  local  coordinate  system,  can  be 
represented  with  respect  to  the  inertial  coordinate  system 
through  the  transformation: 

a„=R:-d„  (8) 

If  we  assume  that  the  receiving  elements  are  constructed  in 
the  same  way  as  the  transmitters,  then  we  can  model  the  722-th 

receive  element,  located  at  position  Xhm ,  as  measuring  the 
received  field  according  to  the  unit-norm  (moment)  direction 
bm  ,  referenced  to  the  (x,y,z)  coordinate  system. 
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The  field  measured  by  the  m- th  receiver  from  the  n- th 
transmitter  was  derived  in  [1]  and  [3].  In  this  paper,  we 
extend  this  forward  model  by  assuming  that  each 
measurement  i  of  the  electric  field  is  function  of  the  receiver 

position  rm ,  receiver  moment  direction  bm ,  and  a 
combination  £:jrwa,  V^eF^J  of  Ns  transmitters 

simultaneously  activated,  represented  with  the  index  set  . 

The  /- th  field  measurement  is  related  with  the  contrast 
function  according  to  the  following  equation: 


bfZ2«G(r-’r»)'a. 

n^Ts 

represents  the  coupling  between  a  combination  s  of 
transmitters  and  the  m- th  receiver.  This  can  also  be  expressed 
using  a  vector  P  (see  [1]  for  details).  Therefore,  the  forward 
model  can  be  approximated  in  discrete  form  by  [  1  ] : 

£  =  +  F  (14) 


=  QnG(rbm,r:)-an+H(i)  +  Z  + 

neTs 

kl\\\v{v')hTm-G(rbm,Y')-  JjQn G(r’,r;)-a( 

D  |_"gIT 

k0  =  0)^IjU(i£(i  (10) 

*  =  (r»’b  »’J)  (U) 

where  Qn  =  jCOIU()Ln  for  a  dipole  of  length  Ln ,  and 

Qn  =  —jCO/Ll{)An  for  a  loop  of  area  An  ,  H  represents  the 

nonlinear  term  (i.e.,  multiple  scattering)  and  Z  represents  the 
noise  term.  Matrix  G  represents  the  dyadic  Green’s  function 
of  the  problem.  Expressions  for  free  space,  half-space, 
layered  Green’s  Dyad  are  presented  in  [8]-[10].  Equation  (9) 
represents  the  forward  model  for  RF  tomography  using  a 
combination  s  of  transmitters  activated  at  the  same  time. 
Both  H  and  Z  are  assumed  unknown,  and  represent  undesired 
perturbations  to  the  received  scattered  field.  Since  tunnels  are 
static  objects  and  the  probing  wavelength  is  large  compared 
to  the  targets,  we  can  reasonably  neglect  these  terms,  and  (9) 
becomes  a  linear  operator  on  V  . 

The  next  logical  step  is  to  collect  other  measurements  by 
varying  Yb ,  bm ,  S  to  form  a  measurement  set: 


In  principle  the  estimation  of  V  can  be  performed  as  follows: 

V  =  L ~'(E-P)  (15) 


where  L  1  can  be  computed  using  Tikhonov  Regularization, 
Truncated  SVD,  Backpropagation,  etc.  (see  [1],  [2]). 

However,  in  real  cases,  one  should  expect  some  uncertainties 

in  G  ,  r“ ,  rbm  ,  bm  and  ;  as  such,  once  can  only  estimate 

the  coupling  effect  via  the  vector  P  and  form  the  image 
using  the  following  inversion  procedure: 

f  =  L (16) 

Unfortunately,  each  i- th  entry  of  P  can  be  up  to  50-60  dB 
higher  in  value  than  the  scattering  terms  in  the  i- th  entry  of 
E  .  As  such,  slight  deviations  P  —  P  can  be  still  higher  than 
the  scattering  signals,  and  further  mitigation  is  required. 

A  way  to  overcome  this  problem  is  to  dramatically  reduce  the 

magnitude  of  P  in  order  to  have  the  deviation  vector  P  —  P 

comparable  with  the  magnitude  of  ES  .  The  reduction  of  P 
can  be  achieved  by  a  proper  selection  of  the  amplitude  and 
phase  of  each  radiating  element,  which  is  the  subject  of  next 
Section. 


£  =  {£(0}  (12> 

Moreover,  the  domain  D  can  be  discretized  into  k  voxels, 
each  one  located  at  r  \  .  Therefore,  we  can  define  the 

contrast  vector  F  =  {F(rf^)J,  where  Vk  corresponds  to  the 

contrast  function  evaluated  at  r  \  .  The  operation  of 
integration  in  (9)  repeated  for  each  i- th  measurement  can  be 
discretized  in  a  matrix  L  operating  on  V  . 

Note  that  the  first  term  in  (9): 


IV.  System  Design 

The  method  we  propose  is  to  design  the  transmit  currents  in 
each  dipole  of  each  antenna  such  that  the  direct-path  signals 
are  heavily  reduced  at  the  receiver  side.  This  design  process 
must  be  repeated  for  each  /-th  measurement.  In  order  to 
provide  robustness  to  uncertainties,  we  seek  current 
distributions  that  produce  low  coupling  over  a  small  area 
surrounding  the  m- th  receiver.  This  task  is  accomplished  by 

considering  a  suitable  set  of  A  points  in  the  neighbor  of  Yb  , 

each  one  positioned  at  a  (small)  distance  Ybs  from  Yb  (see 
Fig.  2).  Furthermore,  we  can  allocate  a  relative  importance  at 
each  point  by  assigning  each  point  with  a  weight  W§  . 
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The  power  of  the  coupling  effect  at  the  receiver  position  Xhm 
and  moment  direction  bm ,  and  in  its  surrounding  A  points 
Xb§  can  be  quantified  via  the  functional: 

2 

•'(».)  =  !>.,  Zb-'e,G(r‘,+r;,r;)  a„  ,(17) 

S= 1  neT  s  2 

By  expanding  (17)  and  using  (8)  we  obtain: 

^K)=  X  (18) 

nJeTs 

G„J  =  t,wsQ.G(<+rs<l)  <19> 

S= 1 

*m=(K<)  (20) 

If  the  coupling  needs  to  be  minimized  independently  from  the 
moment  direction  bm  ,  then  Bw  =  I3x3 . 

Clearly,  the  functional  in  (18)  can  be  recast  in  matrix  form: 

J(dn)  =  dH  K  d  (21) 

Where  the  positive  semidefmite  matrix  K  is  defined  by 

1  'Jj  oU  |  lx  |  ***  iVj  'J  i 

K  =  :  ^  0  (22) 

and  d  represents  a  vector  comprised  of  the  design  currents  in 
the  local  coordinate  systems. 

d  =  [d,  •••  dNf  (23) 

The  optimal  choice  of  the  transmitters  currents  dopt  is  the 

vector  that  minimizes  J .  To  avoid  the  trivial  solution,  we 
constrain  the  total  transmitted  power  to  be  some  positive 
scalar.  Since  dopt  is  invariant  with  respect  to  the  choice  of 
the  power,  we  choose  the  unity  for  convenience. 


<U=argmin[./(dn)] 

(24) 

subject  to:  d^-d^l 

(25) 

The  minimization  of  (21)  using  the  constraint  (25)  can  be 
obtained  analytically  using  the  Lagrangian  multiplier  method. 
The  Lagrangian  is  given  by 

A  =  d//-K-d-2(d// d-l)  (26) 

Imposing  V  A  =  0 ,  we  obtain: 

Kd  =  M  (27) 

Thus,  dopt  is  any  eigenvector  corresponding  to  the  smallest 
eigenvalue  of  K  . 

However,  we  also  wish  to  maximize  the  signal-to-noise 
ration  of  the  scattered  signals,  and  this  will  generally  occur 
only  when  each  transmitter  is  radiating  at  its  maximum 
power.  As  such,  we  prefer  to  solve  the  constrained 
optimization  problem  given  by 

^oPt  -  arg  min  [./  (dB )]  (28) 

subject  to  :  ||d„ \f2  =  df  dn  =  1  Vw  (29) 

This  power  constraint  represents  a  non-convex  set  of 
constraints,  and  the  minimization  can  be  performed  using 
standard  techniques  such  as  sequential  quadratic 
programming  (SQP),  or  interior  point  methods  (IP). 

In  practical  cases,  transmitters  may  be  constructed  using  only 
two  co-located  orthogonal  dipole  (i.e.,  crossed  dipoles),  or 
even  a  single  dipole.  There  are  number  of  ways  to  handle 
these  conditions.  In  this  work,  we  simply  impose  the 
condition 

d€n  =  0  (30) 

when  the  n- th  transmitter  is  a  crossed  dipole,  and 

d,  -  0,  —  0  (31) 

when  the  n- th  transmitter  is  a  single  dipole. 


V.  Simulations 

In  order  to  demonstrate  the  validity  of  our  approach,  we 
simulated  the  imaging  of  belowground  void  structures 
representing  an  underground  facility  using  close-in  sensing. 

The  parameter  of  interest  are  f  —  5  MHz,  £D—9,  and 

(JD  —  5x10  4  [S/m].  The  underground  facility  is 

simulated  using  two  cylinders  of  radius  lm  places  10m 
beneath  the  ground  (see  Fig.  4). 

The  sensors  were  modeled  as  six  transceivers  (i.e.,  co¬ 
located  transmitter  and  receiver)  encircling  the  investigation 
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domain  D  (see  Fig.  4).  All  transmitters  are  considered  to  be 
crossed  dipoles,  each  one  radiating  at  its  maximum  allowed 
power.  The  first  dipole  of  each  transceiver  is  oriented  along 

X ,  while  the  second  dipole  is  oriented  along  y ,  so  that 
simplifies. 

We  computed  the  power  set  of  the  set  of  all  transmitters 
and  for  each  subset  (i.e.,  combination  s  of  transmitters)  we 
optimized  design  currents  for  each  receiver.  Among  all 
possible  combinations,  we  selected  only  10  measurements 
that  provide  the  best  coupling  suppression  at  a  particular 

receiver,  and  at  4  equally  distant  points  (|r^|  =  0.5  [m] ) 

from  the  receiver. 

The  optimization  was  performed  using  the  interior  point 
algorithm  [4]  in  the  MATLAB  optimization  toolbox.  The 
algorithm  was  initialized  with  a  feasible  point,  and  converged 
in  tenths  of  seconds  for  each  design.  The  optimization 
process  was  able  to  produce  current  designs  that  put  nulls  in 
excess  of  50dB  at  the  receiver  locations  (one  such  design  is 
shown  in  Fig.  5). 

For  each  combination  s  and  Rx  pair,  the  optimized 
currents  were  inserted  into  the  FDTD  simulator  GPRMAX 
[6]  and  the  forward  scattered  field  at  the  surface  is  calculated. 
The  field  is  then  collected  from  a  dual  polarized  receiver.  In 
total,  120  samples  are  collected  and  inserted  into  (14).  The 
corresponding  underdetermined  matrix  equation  is  solved 
using  Truncated  SVD  [1],  and  the  contrast  function  was 
estimated.  The  resulting  image  is  shown  in  Fig.  6. 
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Figure  4:  Simulation:  Tx  and  Rx  reside  in  the  same  region. 
Two  tunnels  (in  green)  are  located  10m  belowground. 


Figure  5:  Snapshot  of  the  incident  electric  field  at  the 
surface.  The  receiver  located  at  (-60,  0,  0)  is  experiencing  a 
coupling  below  50dB  from  the  peak  value  in  the  domain  D. 
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Figure  6:  Reconstructed  image  using  120  samples  and 
Truncated  SVD.  Depth  slice  of  10m. 

VI.  Conclusions 

Figure  7  shows  the  reconstructed  image  of  the 
underground  scene.  Using  interferometric  RF  tomography, 
samples  were  collected  without  the  need  of  knowing  the 
background  incident  field  at  the  receivers,  since  the  coupling 
was  mitigated  a  priori.  Clearly,  this  strategy  represents  an 
important  step  toward  a  practical  implementation  of  RF 
tomography  for  real-world  belowground  imaging. 
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